
#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"
#define index_h(i,j) ((i) * N + (j))

void sortArray(float arr[], int ord[], int N) {
    // 初始化ord数组为0到N-1的序列
    for (int i = 0; i < N; i++) {
        ord[i] = i;
    }

    // 使用冒泡排序算法对数组进行排序
    for (int i = 0; i < N - 1; i++) {
        for (int j = 0; j < N - i - 1; j++) {
            if (arr[ord[j]] > arr[ord[j+1]]) { //right
                // 交换ord[j]和ord[j+1]的值
                int temp = ord[j];
                ord[j] = ord[j+1];
                ord[j+1] = temp;
            }
        }
    }
}


/* Auxiliary routines prototypes */
void print_eigenvalues( char* desc, MKL_INT n, float* wr, float* wi , float *wl, int tid);
void print_eigenvectors( char* desc, MKL_INT n, float* wi, float* v,
      MKL_INT ldv );

extern "C" void comp_eig(float *a, int N, float *wr,float *wi, float *vr, int ord[], int tid) {
        /* Locals */
        MKL_INT n = N, lda = N, ldvl = N, ldvr = N, info;
        /* Local arrays */
        //float vl[1];
        float wl[N];
//        printf( "LAPACKE_dgeev (row-major, high-level) Example Program Results\n" );
    int is,js; 
float tmp_sum=0.0;
/*
 FILE *file = fopen("output.txt", "w");
    if (file == NULL) {
        printf("Failed to open file for writing.");
        return;
    } 
    //if(tid==0){ 
      tmp_sum=0.0;
            fprintf(file,"\na=\n");
        for(is=0;is<N;is++){
            for(js=0;js<N;js++){
                fprintf(file,"%e ",a[index_h(is,js)]);
             }
            fprintf(file,"\n");
        }
        //printf("tid=%d,a_sum=%e\n",tid,tmp_sum);
    //}

     fclose(file);
   */ 
        /* Solve eigenproblem */
        info = LAPACKE_sgeev( LAPACK_ROW_MAJOR, 'N', 'V', n, a, lda, wr, wi,
                        NULL, ldvl, vr, ldvr );
        /* Check for convergence */
        if( info > 0 ) {
                printf( "The algorithm failed to compute eigenvalues.\n" );
                exit( 1 );
        }
        /* Print eigenvalues */
        print_eigenvalues( "Eigenvalues", n, wr, wi, wl , tid);
        /* Print left eigenvectors */
        //print_eigenvectors( "Left eigenvectors", n, wi, vl, ldvl );
        /* Print right eigenvectors */
        //print_eigenvectors( "Right eigenvectors", n, wi, vr, ldvr );
        //int N = sizeof(arr) / sizeof(arr[0]); // 数组长度
/* if(tid==12||tid==11){
    printf("tid=%d,wl=\n ",tid);
        for (int i = 0; i < n; i++) {
            printf("%.2e ",wl[i]);
        }
 }*/
        // Sort eigenvalues in ascending order
        sortArray(wl, ord, n); // 调用排序函数
        /*
 if(tid==12){
    printf("wl=\n ");
        for (int i = 0; i < n; i++) {
            printf("%.2e ",wl[i]);
        }
        printf("排序后的数组索引：");
        for (int i = 0; i < n; i++) {
            printf("%d ", ord[i]);
            //printf("%.2e ",wl[ord[i]]);
        }
} */
/*
    printf("lambda=\n ");
    for (int i = 0; i < N; i++) {
        int tmp=0;
        for (int j = 0; j < N; j++) {
            if(wl[j]>wl[i]) tmp++;
        }
        ord[i]=tmp;
        printf("%.2e ",wl[i]);
    }*/
        //exit( 0 );
} /* End of LAPACKE_dgeev Example */

/* Auxiliary routine: printing eigenvalues */
void print_eigenvalues( char* desc, MKL_INT n, float* wr, float* wi , float wl[], int tid) {
        MKL_INT j;
//if(tid==12)        printf( "\n %s\n", desc );
   for( j = 0; j < n; j++ ) {
      if( wi[j] == (float)0.0 ) {
          wl[j]=fabs(wr[j]);
//if(tid==12)         printf( " %6.2f", wr[j] );
      } else {
          wl[j]=sqrt(wr[j]*wr[j]+ wi[j]*wi[j]);
  //      if(tid==12) printf( " (%6.2f,%6.2f)", wr[j], wi[j] );
      }
   }
//   printf( "\n" );
}

/* Auxiliary routine: printing eigenvectors */
void print_eigenvectors( char* desc, MKL_INT n, float* wi, float* v, MKL_INT ldv ) {
        MKL_INT i, j;
//        printf( "\n %s\n", desc );
   for( i = 0; i < n; i++ ) {
      j = 0;
      while( j < n ) {
         if( wi[j] == (float)0.0 ) {
//            printf( " %6.2f", v[i*ldv+j] );
            j++;
         } else {
//            printf( " (%6.2f,%6.2f)", v[i*ldv+j], v[i*ldv+(j+1)] );
//            printf( " (%6.2f,%6.2f)", v[i*ldv+j], -v[i*ldv+(j+1)] );
            j += 2;
         }
      }
//      printf( "\n" );
   }
}
